Import libraries¶
In [1]:
from glob import glob
from natsort import natsorted, ns
import geopy
import folium
from geopy.geocoders import Nominatim
from geopy.extra.rate_limiter import RateLimiter
import rasterio
from rasterio.plot import show
from pyproj import Transformer
from rasterio.windows import Window
import rioxarray
import imageio
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline
Search all tifs and sort files¶
In [2]:
def tif_file(path):
tif_file =[]
# using glob library to get all the file with .tif
files = glob(path,recursive = True)
for file in files:
tif_file.append(file)
# sort files with number in the file (ascending order) using natsorted library
tif_file = natsorted(tif_file, alg=ns.IGNORECASE)
return tif_file
In [3]:
# pwd '/Users/adamtky/Desktop/BeCode/3. Project/Project 3 - 3D Houses'
In [4]:
path = './DSM/**/*.tif'
DSM_tif = tif_file(path)
DSM_tif
Out[4]:
In [5]:
path = './DTM/**/*.tif'
DTM_tif = tif_file(path)
DTM_tif
Out[5]:
Geocoding - single address¶
In [6]:
address = 'Abbaye de Bonne-Espérance'
In [7]:
# enter an address in Belgium
# address = input("Enter an address in Walloon Brabant or Hainaut , Belgium:")
def lat_long(address):
# to get the longtitude and latitude of the address entered
locator = Nominatim(user_agent="myGeocoder")
location = locator.geocode(address)
house_lat_long = [location.latitude, location.longitude]
return house_lat_long
def house_map(func):
# to plot address on a map
house_map = folium.Map(location=func,zoom_start=150)
folium.Marker(location=func,popup=lat_long(address)).add_to(house_map)
house_map.save("house_map.html")
return house_map
house_map(lat_long(address))
# Example address:
# Abbaye de Bonne-Espérance
Out[7]:
Convert longtitude & latitude¶
In [8]:
# longtitude and latitude of the address entered
lat,lon = lat_long(address)
lat,lon
Out[8]:
In [9]:
def transform(lon, lat):
# transform to Belgium 'EPSG:31370' coordinate
transformer = Transformer.from_crs("EPSG:4326", crs_to = 'EPSG:31370' ,always_xy=True)
xx, yy = transformer.transform(lon, lat)
print(f'longtitude:{xx} latitude:{yy}' )
return xx,yy
In [10]:
# coordinate EPSG:31370
xx,yy = transform(lon, lat)
Locate tif for given longtitude & latitude¶
In [11]:
# create all bounding box from tifs
def bounds(tifs):
bounds = []
for i in tifs:
src = rasterio.open(i)
bounds.append(src.bounds)
return bounds
src_bounds = bounds(DSM_tif)
In [12]:
src_bounds[:3]
Out[12]:
In [13]:
scr_path = []
if (xx >= src_bounds[0][0] and xx <= src_bounds[0][2]) & (yy >= src_bounds[0][1] and yy <= src_bounds[0][3]):
scr_path.append(DSM_tif[0])
print(f'The house is in this tif : {DSM_tif[0]}')
elif (xx >= src_bounds[1][0] and xx <= src_bounds[1][2]) & (yy >= src_bounds[1][1] and yy <= src_bounds[1][3]):
print(f'The house is in this tif : {DSM_tif[1]}')
else:
print('The house is not found in the tif')
In [14]:
# check which tif contains the location in DSM
def check_DSM(xx,yy):
scr_path = []
if (xx >= src_bounds[0][0] and xx <= src_bounds[0][2]) & (yy >= src_bounds[0][1] and yy <= src_bounds[0][3]):
scr_path.append(DSM_tif[0])
print(f'The house is in this tif : {DSM_tif[0]}')
elif (xx >= src_bounds[1][0] and xx <= src_bounds[1][2]) & (yy >= src_bounds[1][1] and yy <= src_bounds[1][3]):
print(f'The house is in this tif : {DSM_tif[1]}')
scr_path.append(DSM_tif[1])
else:
print('The house is not found in the tif')
return scr_path
dsm_path = check_DSM(xx,yy)
In [15]:
# check which tif contains the location in DTM
def check_DTM(xx,yy):
scr_path = []
if (xx >= src_bounds[0][0] and xx <= src_bounds[0][2]) & (yy >= src_bounds[0][1] and yy <= src_bounds[0][3]):
scr_path.append(DTM_tif[0])
print(f'The house is in this tif : {DTM_tif[0]}')
elif (xx >= src_bounds[1][0] and xx <= src_bounds[1][2]) & (yy >= src_bounds[1][1] and yy <= src_bounds[1][3]):
print(f'The house is in this tif : {DTM_tif[1]}')
scr_path.append(DTM_tif[1])
else:
print('The house is not found in the tif')
return scr_path
dtm_path = check_DTM(xx,yy)
Clip the location of the house from tif¶
In [16]:
dsm_path
Out[16]:
In [17]:
dtm_path
Out[17]:
In [18]:
def clip_dsm(path,window_size:int):
xds = rioxarray.open_rasterio(path,masked=True,chunks=True)
# set window size
n = window_size
# create coordinates and geometries
coor1,coor2 = [(xx-n),(yy+n)],[(xx+n),(yy+n)]
coor3,coor4 = [(xx+n),(yy-n)] ,[(xx-n),(yy-n)]
geometries = [ {'type': 'Polygon', 'coordinates': [[coor1,coor2, coor3,coor4,coor1 ]]}]
# clip the image as per the geometries size
clipped = xds.rio.clip(geometries)
clip_dsm = clipped.rio.to_raster("clipped_dsm.tif",tiled=True, dtype="int32")
return clipped.plot()
In [19]:
def clip_dtm(path,window_size:int):
xds = rioxarray.open_rasterio(path,masked=True,chunks=True)
# set window size
n = window_size
# create coordinates and geometries
coor1,coor2 = [(xx-n),(yy+n)],[(xx+n),(yy+n)]
coor3,coor4 = [(xx+n),(yy-n)] ,[(xx-n),(yy-n)]
geometries = [ {'type': 'Polygon', 'coordinates': [[coor1,coor2, coor3,coor4,coor1 ]]}]
# clip the image as per the geometries size
clipped = xds.rio.clip(geometries)
clip_dtm = clipped.rio.to_raster("clipped_dtm.tif",tiled=True, dtype="int32")
return clipped.plot()
In [20]:
dsm_path[0]
Out[20]:
In [21]:
dtm_path[0]
Out[21]:
In [22]:
clip_dsm(dsm_path[0],50)
Out[22]:
In [23]:
clip_dtm(dtm_path[0],50)
Out[23]:
In [24]:
def chm_tif():
# open the digital terrain model
with rasterio.open('clipped_dtm.tif') as src:
lidar_dtm_im = src.read(1, masked=True)
dtm_meta = src.profile
# open the digital surface model
with rasterio.open('clipped_dsm.tif') as src:
lidar_dsm_im = src.read(1, masked=True)
dsm_meta = src.profile
# calculate canopy height model
lidar_chm = lidar_dsm_im - lidar_dtm_im
with rasterio.open('clipped_chm.tif', 'w', **dsm_meta) as ff:
ff.write(lidar_chm,1)
chm_tif = 'clipped_chm.tif'
return chm_tif
Create 3D¶
In [25]:
def House_3D(tif,size:int,azim=215):
chm = imageio.imread(tif)
nx,ny = chm.shape
x = np.linspace(0, size*2, nx)
y = np.linspace(0, size*2, ny)
yv,xv = np.meshgrid(x, y)
fig = plt.figure(figsize=(15,15))
ax = fig.add_subplot(111, projection='3d')
chm3d=ax.plot_surface(xv,yv,chm,cmap='plasma',linewidth=0)
ax.set_title('CHM(Canopy Height Model)')
ax.set_xlabel('Distance (m)')
ax.set_ylabel('Distance (m)')
ax.set_zlabel('Elevation (m)')
ax.view_init(azim=azim)
fig.colorbar(chm3d, shrink=0.3, aspect=10);
fig.savefig('house.png', dpi=200)
return plt.show()
In [26]:
House_3D(chm_tif(),50,300)
In [27]:
# compare google map image of the house
from IPython.display import Image
Image("Basilique Notre-Dame de Bonne-Espérance.png")
Out[27]:
In [28]:
# to view actual house on google map
import webbrowser
url = 'https://www.google.com.my/maps/place/'+str(lat)+','+str(lon)
webbrowser.open(url)
Out[28]:
In [29]:
from osgeo import gdal
from mayavi import mlab
ds = gdal.Open(chm_tif())
data = ds.ReadAsArray()
data = data.astype(np.float32)
mlab.figure(size=(640, 800), bgcolor=(0.16, 0.28, 0.46))
surf = mlab.surf(data,warp_scale='auto')
#mlab.show()
In [30]:
Image("Basilique Notre-Dame de Bonne-Espérance_mayavi.png")
Out[30]:
